import struct
import math
import numpy as np
from matplotlib import pyplot

def anint(x, lx):
    lhx=0.5*lx;
    while x>lhx:
        x-=lx
    while x<-lhx:
        x+=lx
    return x

def distancia(lbox,pos1,pos2):
    dis2=0
    drx=pos2[0]-pos1[0]
    drx=anint(drx,lbox[0])
    dry=pos2[1]-pos1[1]
    dry=anint(dry,lbox[1])
    drz=pos2[2]-pos1[2]
    drz=anint(drz,lbox[2])
    dis2=drx*drx+dry*dry+drz*drz
    return  [math.sqrt(dis2),drx,dry,drz]

def Lboundary(x, lx):
    res=x
    while x>=lx:
        x-=lx
    while x<0:
        x+=lx
    return x

def boundary(pos, lbox):
    pos[0]=Lboundary(pos[0],lbox[0])
    pos[1]=Lboundary(pos[1],lbox[1])
    #pos[2]=Lboundary(pos[2],lbox[2])
    return pos

def binaryExtract(f):
    s = f.read(8)
    t= struct.unpack("d", s)[0]
    s = f.read(4)
    nacc = struct.unpack("i", s)[0]
    s = f.read(4)
    npart = struct.unpack("i", s)[0]
    s = f.read(8*3)
    lbox= struct.unpack("ddd", s)
    posiciones=[]
    velocidades=[]
    fuerzas=[]
    diametros=[]
    #Comienza el bucle sobre el numero de particulas
    for i in range(npart):
        #Extraemos las posiciones
        s = f.read(8*3)
        pos= list(struct.unpack("ddd", s))
        #Extraemos las velocidades
        s = f.read(8*3)
        vel= struct.unpack("ddd", s)
        #Extraemos las fuerzas
        s = f.read(8*3)
        force= struct.unpack("ddd", s)
        #Extraemos el diametro
        s = f.read(8)
        diam= struct.unpack("d", s)[0]
        boundary(pos,lbox)
        posiciones.append(pos)
        velocidades.append(vel)
        fuerzas.append(force)
        diametros.append(diam)
    return [npart, nacc, lbox, posiciones,velocidades,fuerzas,diametros]



def conectividad(pos, lbox, rcontact, zseg):
    npart = len(pos)
    angulos=[]
    #Recorremos todos los pares de particulas
    for i in range(npart-1):
        for j in range(i+1,npart):
            [dis,dx,dy,dz] = distancia(lbox,pos[i],pos[j])
            if dis < rcontact and pos[i][2]<(lbox[2]-zseg) and pos[i][2]>zseg and pos[j][2]<(lbox[2]-zseg) and pos[j][2]>zseg:
                #Las particulas estan en cntacto-> extraemos theta entre i-j
                theta=math.acos(dz/dis)*180.0/math.pi
                angulos.append(theta)
                #Por simetria incluimos angulo entre j-i
                angulos.append(180.0-theta)
    return angulos

def histogram(minValue,maxValue,values, nbins):
    suma=0.0
    toRad=math.pi/180.0;
    dist = np.zeros((nbins))
    count = np.zeros((nbins))
    step=(maxValue-minValue)/float(nbins)
    print step,nbins
    for val in values:
        if val>=minValue and val<=maxValue:
            dnd=int(math.floor((val-minValue)/step))
            dist[dnd]+=val
            count[dnd]+=1.0
            suma+=1.0
    
    x=[]
    y=[]
    sss=0.0
    for i in range(len(count)):
        x.append(minValue+float(i+0.5)*step)
        fracArea=(math.cos(toRad*(x[i]-0.5*step))-(math.cos(toRad*(x[i]+0.5*step))))*0.5
        y.append(count[i]/fracArea/suma)
    
    #Normalizacion
#    suma=0.0
#    for i in range(1,len(y)):
#        suma+=step*(y[i]+y[i-1])*0.5

        
    return [x,y]


rcontact = 1.05
nbins=180
zseg=2.0
name="/home/jcfernandez/remote/rperi76/simulaciones/eps1T0.1/u100/AMPL/W2.0/A0.05/mini_config.bin"
f = open(name,'r')
[npart, nacc, lbox, posiciones,velocidades,fuerzas,diametros]=binaryExtract(f)
angulos = conectividad(posiciones, lbox, rcontact,0.0)
[x,y]=histogram(0.,180.,angulos, nbins)
pyplot.plot(x,y,'ro-')
angulos = conectividad(posiciones, lbox, rcontact,1.0)
[x,y]=histogram(0.,180.,angulos, nbins)
pyplot.plot(x,y,'go-')
angulos = conectividad(posiciones, lbox, rcontact,3.0)
[x,y]=histogram(0.,180.,angulos, nbins)
pyplot.plot(x,y,'bo-')
pyplot.show()
f.close()